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Abstract 

We propose an algorithm based on Newton's method and subdivi- 
sion for finding all zeros of a polynomial system in a bounded region 
of the plane. This algorithm can be used to find the intersections 
between a line and a surface, which has applications in graphics and 
computer-aided geometric design. The algorithm can operate on poly- 
nomials represented in any basis that satisfies a few conditions. The 
power basis, the Bernstein basis, and the first-kind Chebyshev basis 
are among those compatible with the algorithm. The main novelty 
of our algorithm is an analysis showing that its running is bounded 
only in terms of the condition number of the polynomial's zeros and 
a constant depending on the polynomial basis. 

1 Introduction 

The problem of line-surface intersection has many applications in areas such 
as geometric modeling, robotics, collision avoidance, manufacturing simula- 
tion, scientific visualization, and computer graphics. It is also a basis for con- 
sidering intersections between more complicated objects. This article deals 
with intersections of a line and a parametric surface. The parametric method 
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of surface representation is a very convenient way of approximating and de- 
signing curved surfaces, and computation using parametric representation is 
often much more efficient than other types of surface representations. 

Typically, intersection problems reduce to solving systems of nonlinear 
equations. Subdivision methods introduced by Whitted [221 US] were the 
first to be used for this problem. In these methods, a simple shape such as 
rectangular box or sphere is used to bound the surface and is tested whether 
the line intersects the bounding volume. If it does, the surface patch is 
subdivided, and the bounding volumes are found for each subpatch. The 
process repeats until no bounding volumes intersect the line or the volumes 
are smaller than a specified size where the intersection between such volumes 
and the line are taken as the intersections between the surface and the line. 
Subdivision methods are robust and simple, but normally not efficient when 
high accuracy of the solutions are required. They also cannot indicate if 
there are more than one zero inside the remaining subdomains. 

Regardless, a variation of subdivision methods known as Bezier clipping 
by Nishita et al. should be noted for its efficient subdivision [13J. For a 
non-rational Bezier surface, Bezier clipping uses the intersection between 
the convex hull of the orthographic projection of the surface along the line 
and a parameter axis to determine the regions which do not contain any 
intersections before subdividing the remaining region. Sherbrooke and Pa- 
trikalakis generalizes Bezier clipping to a zero-finding algorithm for an n- 
dimensional nonlinear polynomial system called Projected Polyhedron (PP) 
algorithm [T6] . 

On the numerical side, Kajiya [TO] proposes a method for intersecting a 
line with a bicubic surface based on algebraic geometry without subdivisions. 
His method is robust and simple. However, it is too costly to extend to higher 
degree polynomials. Jonsson and Vavasis [H] introduce an algorithm for solv- 
ing systems of two polynomials in two variables using Macaulay resultant 
matrices. They also analyze the accuracy of the computed zeros in term of 
the conditioning of the problem. 

Another approach is to combine a subdivision method with a Newton- 
type method, using the latter to find the solutions of the resulted system of 
equations after subpatches pass some criteria. The advantages of Newton's 
method are its quadratic convergence and simplicity in implementation, but 
it requires a good initial approximation to converge and does not guarantee 
that all zeros have been found. To remedy these shortcomings, Toth [20] 
uses the result from interval analysis to determine the "safe regions" , where 
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a Newton-like method starting from any point in them are guaranteed to 
converge. He tests each patch if it is a safe region and if its axis-aligned 
bounding box intersects the line. If neither is true, the patch is subdivided 
recursively. Lischinski and Gonczarowski [12] propose an improvement to 
Toth's algorithm specific to scene-rendering in computer graphics by utilizing 
ray and surface coherences. 

In contrast, other researchers develop methods to estimate good initial 
points for Newton's method rather than test for convergence of each choice 
of initial points. These methods use tighter bounding volumes and subdivide 
the surface adaptively until subpatches are flat enough, that is, until they 
are close enough to the bounding volumes. Then the intersection between 
the bounding volume and the line is chosen as the initial point for Newton's 
method. Examples of methods in this category are [TJ [6j [HI [181 123] . There 
is also the ray-tracing algorithm by Wang et al. that combines Newton's 
method and Bezier clipping together [2T| . 

Our algorithm is in the same category as Toth's in that it tests for the 
convergence of initial points before performing Newton's iteration to find 
solutions and it uses a bounding polygon of a subpatch to exclude the one 
that cannot have a solution. The convergence test our algorithm uses is 
the Kantorovich test. Because the Kantorovich test also tells us whether 
Newton's method converges quadratically for the initial point in question in 
addition to whether it converges at all, we can choose to hold off Newton's 
method until quadratic convergence is assured. 

The main feature of our algorithm is that there is an upper bound on 
the number of subdivisions performed during the course of the algorithm 
that depends only on the condition number of the problem instance. For 
example, having a solution located exactly on the border of a subpatch does 
not adversely affect its efficiency. To the best of our knowledge, there is 
no previous algorithm in this class whose running time has been bounded 
in terms of the condition of the underlying problem instance, and we are 
not sure whether such an analysis is possible for previous algorithms. Its 
efficiency also depends on the choice of basis because the type of bounding 
polygon depends on the basis. 

The notion of bounding the running time of an iterative method in terms 
of the condition number of the instance is an old one, with the most notable 
example being the condition- number bound of conjugate gradient (see Chap- 
ter 10 of (8]). This approach has also been used in interior-point methods for 
linear programming [7] and Krylov-space eigenvalue computation [T9] . 
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2 The theorem of Kantorovich 



Denote the closed ball centered at x with radius r > by 

B(x, r) = {y e R n : \\y - x\\ < r}, 

and denote B(x, r) as the interior of B(x, r). Kantorovich's theorem in affine 
invariant form is 

Theorem 2.1 (Kantorovich, affine invariant form 0[TT]). Let f : D C R n — > 

M. n be differentiable in the open convex set D. Assume that for some point 
x° G D, the Jacobian f'(x°) is invertible with 

IIA* )- 1 /^ )!! < v . 

Let there be a Lipschitz constant uo > for f'{x°)~ 1 f' such that 

\\f(x°)-\f'(x) - f'(y)) || < u ■ \\x - y\\ for all x,y E D. 
If h = rjuj < 1/2 and B(x°,p-) C D, where 

1 - VI - 2h 

P- = , 

to 

then f has a zero x* in i?(x°,p_). Moreover, this zero is the unique zero of 
f in {B(x°, p_) U 5(2°, p+)) n D where 



1 + VI - 2h 

P+ = 

u 

and the Newton iterates x k with 

x k+i = x k_ f'^yif^ 

are well-defined, remain in B(x°,p-), and converge to x* . In addition, 

We call x° a fast starting point if the sequence of Newton iterates start- 
ing from it converges to a solution x* and ([1]) is satisfied with h < 1/4, 
which implies quadratic convergence of the iterates starting from x°. The 
Kantorovich's theorem also holds for complex functions [5J. 
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3 Formulation and representation of the line- 
surface intersection problem 

Let 0o, . • • , <p n denote a basis for the set of univariate polynomials of degree 
at most n. For example, the power basis is defined by <pi(t) = t i . Other 
choices of basis are discussed below. 

Let S be a two-dimensional surface embedded in M 3 parametrized by 

f(u,v) = E^oE™=o^j<M«)0j(f), 0<u,v<l, 

where G MP (i = 0, 1, . . . , m; j — 0, 1, . . . , n) denote the coefficients. Define 
a line 

L = {p + dt: t G E}, 

where p, d G M 3 , d ^ 0. The line-surface intersection problem is to find all of 
the intersections between 5" and L, which are the solutions of the polynomial 
system 

f{u,v)-{p + dt) = 0. (2) 

The system fl2]) can be reduced to a system of two equations and two un- 
knowns. To show this, we first break (J2J) into its component parts 

/i(it,f) -pi -tdi = 0, 
f 2 (u,v) -p 2 -td 2 = 0, 
f 3 (u,v) -p 3 -td 3 = 0. 

Here, the subscript i denotes the ith coordinate of the point in three-dimensional 
space. Assuming \di\ > max{|<i 2 |, \d 3 \}, we have the equivalent system 

di(h(u,v) -p 2 ) - d 2 (fi(u,v) -px) = 0, , 3 - 
di(f 3 (u,v) -p 3 ) - d 3 (fi(u,v) -pi) = 0, 

which can be rewritten with the same basis <fii(u)(f)j{v) (see item H] on the list 
of basis properties below) as 

m n 

f(u, v) = J2J2 C >;M U )M V ) = °- ( 4 ) 

i=0 j=0 

The system fll]) is the one our algorithm operates on. 
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Since the parametric domain of the surface under consideration is square, 
our algorithm uses the infinity norm for all of its norm computation. There- 
fore, for the rest of this article, the notation ||-|| is used to refer specifically 
to infinity norm. 

Our algorithm works with any polynomial basis (f>i(u)<f>j(v) provided that 
the following properties hold: 

1. There is a natural interval [I, h] that is the domain for the polynomial. 
In the case of Bernstein polynomials, this is [0, 1], and in the case of 
power and Chebyshev polynomials, this is [—1, 1]. 

2. It is possible to compute a bounding polygon P of S — {f(u,v) : I < 

u,v < h}, where f(u,v) = Y^Lo S"=o c ij<l>i( u ) ( f ) j( v ) and c ij e K<i for 
any d > 1, that satisfies the following properties: 

(a) Determining whether G P can be done efficiently (ideally in 
0{mn) operations). 

(b) Polygon P is affinely and translationally invariant. In other words, 
the bounding polygon of {Af(u, v) + b : I < u, v < h} is {Ax + b : 
x G P} for any nonsingular matrix A G M dxd and any vector 
b G R d . 

(c) For any y G P, 

\\y\\<6 max ||/(«,t;)||, (5) 

l<u,v<h 

where 9 is a function of m and n. 

(d) If d = 1, then the endpoints of P can be computed efficiently 
(ideally in 0(mn) time). 

3. It is possible to reparametrize with [/, h} 2 the surface Si = {f(x) : x G 
B(x°,r)}, where x° G M 2 and r G K. > 0. In other words, it is possible 
(and efficient) to compute the polynomial / represented in the same 
basis such that Si = {/(£) : x G [I, h] 2 }. 

4. Constant polynomials are easy to represent. 

5. Derivatives of polynomials are easy to determine in the same basis, 
(preferably in 0(mn) operations). 
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Recall that P is a bounding polygon of S if and only if x G S implies igP. 

As shown later in Section [3, the efficiency of our algorithm depends on 
0. Hence, the choice of the basis affects the algorithm's performance as each 
basis allows different ways to compute bounding polygons. 

4 The Kantorovich-Test Subdivision algorithm 

Before we detail our algorithm, we define notation and crucial quantities that 
are used by the algorithm and its analysis. Denote x = (u, v) as a point in 
two-dimensional parametric space and f(x) = f(u,v) as the value of / at x. 

For a given zero x* of polynomial /, let uj*(x*) and p*(x*) be quantities 
satisfying the conditions that, first, u>*(x*) is the smallest Lipschitz constant 
for fix*)- 1 /', i-e., 

||f (xT 1 (/'(*) - f'(y))\\ < ■ 11^ — 2/11 for all x,y G B(x*,p*(x*)) 

(6) 

and, second, 

p*(**) = — 7—r- (7) 

Since u^x*) is nondecreasing as p*(x*) increases in (jHJ) but p*(x*) is decreas- 
ing as uJi,{x*) increases in (JTj), there exists a unique pair (u*(x*) , p*(x*)) sat- 
isfying the above conditions, and this pair can be obtained by binary search. 
When more than one function is being discussed, we use uj((x*) to denote 
uj*(x*) of the function /. Approximation to these two quantities, uj*(x*) and 
p*(x*), are computed and made use of by the algorithm. 

For clarity, we simply abbreviate uj^x*) and p*(x*) as and p*, respec- 
tively, throughout the rest of this article when it is clear from the context to 
which x* the quantities belong. 

Straightforward application of the affine invariant form of Kantorovich's 
theorem with x° = x* and D = B( *)) yields the result that x* is the 

unique zero of / in B(x*, p*(x*)). In fact, the above definitions of uj*{x*) and 
p*(x*) are chosen such that the ball that is guaranteed by the Kantorovich's 
theorem to contain no other zeros than x* is the largest possible. 

Define 

7(0) =V (V0(40 + 1) - 80) , 

where 6 is as in ©. Observe that 1 < 7 (0) < (V5 + 2) /4 w 1.0590 since 
7(0) is a decreasing function for positive and > 1 by the definition of a 
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bounding polygon. Another quantity of interest is ojd 1 , which is defined as 
the smallest nonnegative constant u satisfying 



ll/'^T 1 (/'(</) -/'(*)) II < w\\v-z\\, y,z e D',x* e[o,i] 2 

satisfying f(x*) = 0, 



(8) 



where 

/y = [- 7 (£),l + 7 (£)] 2 . (9) 

The motivation of this definition of D' is that it contains all domains whose 
Lipschitz constants may be needed during the course of the algorithm. De- 
note ujf as the maximum of ud> and all uj*(x*) 

ojf = max{u;z)/, max ljJx*)}. 

x*£C 2 :f(x*)=0 

Finally, define the condition number of / to be 

cond( /) = max{u; f , max \\ f (x*)- 1 f (y)\\} . (10) 

Note that in (181), x* is restricted to zeros in [0, l] 2 whereas in (fTOj) . x* ranges 
over all complex zeros of /. We defer the discussion of why (II Op is a reasonable 
condition number until after the description of our algorithm. 

We define the Kantorovich test on a region X = B(x°,r) as the appli- 
cation of Kantorovich's Theorem on the point x° using B(x°,2-f(e)r) as the 
domain D in the statement of the theorem and ||/'(:e ) _1 /(:e )|| as V- F° r 
lu, we instead use Gj > u, where Cj is defined by ([Ti|) below. The region X 
passes the Kantorovich test if t]Cj < 1/4 and B(x°,p^) C D', which implies 
that a; is a fast starting point. 

The other test our algorithm uses is the exclusion test. For a given region 
X, let fx be the polynomial in the basis (pi(u)<j)j(v) that reparametrizes with 
[I, h} 2 the surface defined by / over X. The region X passes the exclusion 
test if the bounding polygon of {fx(u,v) : I <u,v < h} excludes the origin. 
Note that the bounding polygon used in this test must satisfy item [2] of the 
basis properties listed in Section [31 

We now proceed to describe our algorithm, the Kantorovich- Test Subdi- 
vision algorithm or KTS in short. 

Algorithm KTS: 

• Let Q be a queue with [0, l] 2 as its only entry. Set 5 = 0. 
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• Repeat until Q = 

1. Let X be the patch at the front of Q. Remove X from Q. 

2. If X £ X 5 for all X 5 e S, 

— Perform the exclusion test on X = B(x°,r) 

— If X fails the exclusion test, 

(a) Perform the Kantorovich test on X 

(b) If X passes the Kantorovich test, 

i. Perform Newton's method starting from x° to find a 
zero x*. 

ii. If x* £~ X s for any X s £ S (i.e., x* has not been found 
previously), 

* Compute p*(x*) and its associated ou*(x*) by binary 
search. 

* Set S = SU {B(x*,p,(x*))}. 

(c) Subdivide X along both u and w-axes into four equal sub- 
regions. Add these subregions to the end of Q. 

A few remarks are needed regarding the description of the KTS algorithm. 

• The subdivision in step 2.c is performed regardless of the result of the 
Kantorovich test. In general, passing the Kantorovich test does not 
imply that there is only one zero in X. 

• The check that the zero found by Newton's method is not a duplicate 
(step 2.b.ii) is necessary since the Kantorovich test may detect a zero 
outside X. 

• If the Kantorovich test is not applicable for a certain patch due to the 
Jacobian of the midpoint being singular, the patch is treated as if it 
fails the Kantorovich test. 

One property of KTS is that it is affine invariant. In other words, left- 
multiplying / with a 2-by-2 matrix A prior to executing KTS does not change 
its behavior. This is the main reason we define the condition number to be 
affine invariant. Define g = Af. To see that our condition number is affine 
invariant, note that g\x)- l g\y) = [Af\x)Y l Af'(y) = f'(x)~ 1 A~ 1 Af'(y) = 
f'{x)~ 1 f'(y) for any x,y £ W 1 . Therefore, cond(g) = cond(/). In contrast, 
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simpler condition numbers such as Lipschitz constants for /' are not affine 
invariant and hence are not chosen for our analysis. 

Since Toth's algorithm is the most similar one to KTS, it is worthwhile 
to discuss the main differences between the two and the implications these 
differences make. First, Toth's uses the Krawczyk-Moore test and another 
unnamed test, both based on interval analysis, as the convergence test. These 
two tests guarantee linear convergence for the simple Newton iteration — 
a variation of Newton's method where the Jacobian of the initial point is 
used in place of the Jacobian of the current point in every iteration. With 
our Kantorovich test, KTS starts Newton's method only when quadratic 
convergence is assured. 

Another main difference is in the choice of domains for the convergence 
test. Toth's uses the subpatch X itself as the domain for the test. This 
choice may exhibit undesirable behavior when a zero lies on the border of a 
subpatch, which is not necessarily on or near the border of the original domain 
[0, l] 2 . For example, consider the function f(u, v) = (it 2 — .25, v — .8) T whose 
zeros are (.5, .8) and (—.5, .8). The patch {(it, v) : .5 < u < .5 + e, a < v < b} 
does not pass either of Toth's convergence tests for any e > and any 
a < .8 < b although the patch {(it, v) : .45 < u < .8, < v < 1}, a large 
patch whose borders do not coincide with any zeros, does pass the Krawczyk- 
Moore test. This results in excessive subdivisions by Toth's algorithm. KTS 
uses B(x°, 27(0)r) as the domain for X to avoid this problem. Theorem 17.11 
below shows that the Kantorovich test does not have trouble detecting the 
zeros locating on the border of the subpatch. 



5 Implementation details when using power, 
Bernstein, or Chebyshev bases 

This section covers the implementation details of KTS when the polynomial 
system is in the power, Bernstein, or Chebyshev bases. The power basis for 
polynomials of degree n is <fik{t) = t k (0 < k < n). The Bernstein basis is 

(j) k (t) = Z k>n (t) = ^ j (1 - t) n -H k (0 < k <n). The Chebyshev basis is 

4>k{t) = Tkit) (0 < k < n), where T^it) is the Chebyshev polynomial of the 
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first kind generated by the recurrence relation 



T k+1 (t) 



t 



2tT k {t)-T k _ x {t) for k> 1. 



1. 



(11) 



5.1 Bounding polygons 

We begin with the choices of / and h and the definitions of bounding polygons 
of the surface S = {f(u,v) : I <u,v < h}, where f(u,v) is represented by 
one of the three bases, that satisfy the required properties detailed in Section 
|3j For Bernstein basis, the convex hull of the coefficients (control points), 
call it Pi, satisfies the requirements for I = and h = 1. The convex hull Pi 
can be described as 



satisfies the requirements for / = — 1 and h = 1. Note that is a bounding 
polygon of S in the Chebyshev case since |Tfc(£)| < 1 for any k > and any 
t G [—1,1]. Determining whether G P2 is done by solving a small linear 
programming problem. The value of 9 for each case is summarized in Table 
[TJ The reader is referred to [17] for the derivation of 9 for each of the three 
bases as well as the proofs that Pi and P2 satisfy all of the basis properties 
listed in Section [31 

5.2 Computation of Lipschitz constant 

Another step of KTS that needs further elaboration is the computation 
of Lipschitz constant in the Kantorovich test. The Lipschitz constant for 
/'(x ) -1 /' = g is obtained from an upper bound of the derivative of g 




(12) 



For power and Chebyshev bases, the bounding polygon 




(13) 
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Basis 


e 


Bernstein 

Chebyshev 
Power 


( ^mii y~\ max{ \m— %' |, \i'\} 

{2^1=0 Ik'^i \i-i>\ j 
{m + l)(r 


(\^ n TT m^{\n~j'\,\j'\}\ 

\2^j oil,',-., \j-f\ J 
2(m+l)(n+l) 
i + l)(3 m+1 - l)(3 n+1 - 1)/ 


= 0(m m+1 n" +1 ) 

2 



Table 1: The value of #'s of the power, the Bernstein, and the Chebyshev 
bases and their corresponding bounding polygons. 

for all x G X. Let g = gx he the polynomial in the same basis as / that 
reparametrizes with [I, h} 2 the surface defined by g over X. We have that 

maxims) || = max ||#'(x)|| 

xeX x£[l,h] 2 

= max max ||^'(x)y|| 

xe[i,h] 2 II2/IN 1 

2 2 

< max max Yl ' ( x ) ' 

xG \l,h] 2 i z — ' z — ' 
1 J j=l k=l 

< 4 max max \g' ijk (x)\. 

i,j,k xe[l,h] 2 J 

Note that each entry of g' can be written as a polynomial in the same basis 
as / (refer to property of the basis). For this reason, an upper bound of 
max se[i,A] 2 \9ijk( x )\ can be computed as follows: Let P^ be the bounding 
polygon (bounding interval in this case) of {g'ij k (x) : x G computed 
in the same way as described in Section 15. 1[ The maximum absolute value 
of the endpoints of P^ (refer to property [2d] of the basis) is an upper bound 
of max l6 p iA ]2 |g^- fc (a;)|. Let Co denote the Lipschitz constant computed in this 
manner, that is, 

Co = 4 max max | g ijk (x)\, (14) 

i,3,k x&{l,h] J 

where max^g^^p \g[j k {x)\ is computed from the endpoints of its bounding 
interval. 



6 Significance of our condition number 

We now discuss the significance of ffTUl) to the conditioning of the problem. 
In particular, we attempt to justify that the efficiency of any algorithm in 
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the same class as KTS is dependent on fllOp . This class of algorithms being 
considered includes any algorithm that (i) isolates unique zeros with subdi- 
vision before finding them and (ii) will not discard a patch until the convex 
hull of its function values (which is clearly a subset of any possible bounding 
convex polygon) excludes the origin. 



6.1 Condition number and the Kantorovich test 

This section discusses the relationship between uif and the Kantorovich test. 
We show that, for any given zero x* of an arbitrary /, there is a function 
/ such that x* is also a zero of /, f'(x*) = f(x*), uj((x*) = uj((x*), and / 
has another zero y* with \\y* — x*\\ = p*(x*). For example, consider a zero 
x* = (.5, .5) of the function / = (u 3 - 2.2u 2 + l.hhu - .35, v 2 - .7v + .if, 
of which = .1. A corresponding / with the above properties is / = 

{u 2 — .9u + .2,.3v — .15) T , which has zeros at (.5, .5) and (.4, .5). Since 
the Kantorovich test uses only the function value, its first derivative, and the 
Lipschitz constant, all of which are the same for / and / at x*, the functions / 
and / are identical from the perspective of the Kantorovich test applied to x*. 
Therefore, p*(x*) is a reasonable number that quantifies the distance between 
x* and its nearest other zero barring the usage of additional information. 
Consequently, uif, which is greater than or equal to u*(x*) = 2/p*(x*) for 
all zeros x* of /, describes the distance between the closest pairs of zeros 
of /. Therefore, the efficiency of any algorithm that isolates unique zeros is 
dependent on Uf. 

The function / with the above properties can be constructed as follows: 
Let x* = (u*,v*), f'(x*) = ( ai a<1 ), and uj((x*) = uj. If [« 4 | > \a 3 \, 

\ «3 «4 / 



(u — u*) 2 + a>i(u — u*) + a 2 (v — v*) 



f(u,v) = 2a 4 V \,f\ A ' ■ 15 

Otherwise, 

fU v) = ( a ^ u ~ u *) + ^S^{v ~ v*) 2 + a 2 (v - 
\ 013(11 — It*) + oti(v — V*) 

It is straightforward to verify that f(x*) = 0, f'(x*) = f'(x*), and 
uj{(x*) = 00. We now show that \\y* — x*\\ = p*(x*) for the case where 
1^4 1 > I c^3 1 . The other case can be verified in the same manner. Let 
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y* = (u* + Au,v* + Av). Substituting y* into (TTBT) and setting it to zero 
yields 



g(u* + Au,v* + Av) 



Solving (fi~6l) yields 



An = --, 

to 

Av = — • — . 

«4 tO 

Since \a±\ > \a 3 \ and p*(x*) = 2/to, \\y* — x*\\ = p*(x*). 



'^{Auf + ai Au + a 2 Av \ = f \ 
a 3 Au + a 4 Av ) ~ \ )' 

(16) 



6.2 Condition number and the exclusion test 

The other term in our condition number, max 3 .. eC 2._ f ( 3 ,.) =0i!/e [ 0il ]2 H/'Oe*) -1 /' 
relates to the convex bounding polygon test — the test to determine whether 
the convex bounding polygon of a subpatch contains the origin. We show 
that there exists a function / such that a patch B(x°, r) where x° is relatively 
close to a zero, fails the convex bounding polygon test if r > 0(1/ cond(/)). 
Denote x° = (u°, v°). Define the complex function g(z) = (z — (u° — e — ie))- 
(z — (u° + e — ie)), where eel and < e < 1. Consider the following func- 
tion 

f(u,v) = 



Re (g(u + iv)) 
Im (g(u + iv)) 

u 2 -v 2 - 2u°u - 2ev - 2e 2 + (u ) 2 
2uv - 2u°v + 2eu - 2eu° 



(17) 



where Re(^) and lm(z) denote the real and imaginary parts of the complex 
number z, respectively. The four complex zeros of / are (u° — e, — e) , (u° + 
e, — e), (u°, — e — ie), and (u°, —e + ie). Therefore, 

max ||/Vr7'(y)|| = C>(l/e). 

x*eC 2 :/(a;*)=0,yG[0,l] 2 



Moreover 



u>, = 0(l/e). 
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We now show for the case that v° = 0(e) that B(x°,r) fails the convex 
bounding polygon test if r > 0(e). Let A be the circular arc centered at 
(u°, —e) that goes from (u° + r, v — r) to (u° — r,v° — r) counterclockwise. 
Observe that / maps A to the circular arc centered at (— e 2 ,0) that goes 
from (2(v° + e)r - 2ev° - (v ) 2 - 2e 2 , 2r(v° + e - r)) to (2(v° + e)r - 2ev° - 
(v ) 2 — 2e 2 , —2r(v° + e — r)) counterclockwise (see Figure [T]). Notice that 
2r(f° + e — r) > because B(x°,r) C [0, l] 2 . Therefore, the convex bounding 
polygon of f(A) contains the origin if r > ((v ) 2 + 2ev° + 2e 2 )/(2(v° + e)) = 
0(e) (recall the assumption that v° = 0(e)). Since A C B(x°,r), the convex 
bounding polygon of f(B(x°,r)) also contains the origin and the convex 
bounding polygon test fails. 



7 Time complexity analysis 

In this section, we prove a number of theorems relating to the behavior of the 
KTS algorithm. We analyze the efficiency of KTS by showing that a patch 
either is a subset of a safe region, passes the Kantorovich test, or passes 
the exclusion test when it is smaller than a certain size that depends on the 
condition number of the function. Hence, we have the upper bound of the 
total number of patches examined by KTS in order to solve the intersection 
problem. 

Recall that the Lipschitz constant £u given by (|T%|) is not the smallest 
Lipschitz constant of f'(x°)^ 1 f over D', where D' is given by (J9j). However, 
we can show that d) < 48uj, where u denotes the smallest Lipschitz constant 
of /'(x ) -1 / over D' . Since Cj is computed from the endpoints of the bounding 
intervals of max^^^p \g' ijk (x)\, by ©, 

to < AO max max oL t (a;)| 
~ i,j,k xe[i,h] 2 1 3 1 

= 4^maxmax |^ - fc (a;)| 

i,j,k x£X 

< 49 max | \g'(x)\\ = AOu. (18) 

With this bound on &>, we can now analyze the behavior of the Kantorovich 
test. 

Theorem 7.1. Let x° be a point in [0, l] 2 such that f'(x°) is invertihle. 
Let x* be a zero of f that is contained in B(x°,r), where r is the radius of 
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--B(x°,r) 
— A 



vO-r 




uO-r 



uO 
u 



uO+r 



(a) The circular arc A C B(x°,r) . 



0.2p 
0.15 

0.1 
0.05 


-0.05 

-0.1 - 

-0.15 

-0.2 — 
-0.2 



- - convex bounding polygon of f(A) 
— f(A) 




-0.1 



0.1 



0.2 



(b) The range f(A) and its convex bounding polygon 

Figure 1: The circular arc A centered at (u°, — e) that goes from (u°+r, v°—r) 
to (u° — r, v° — r) and its range f(A) where / is as in (IT7|) . Figure [TBI shows 
that the bounding convex polygon of f(A) contains the origin, and therefore 
B(x°,r) fails the convex bounding polygon test. 
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the patch under consideration. The patch X = B(x°,r) C [0, l] 2 passes the 
Kantorovich test if 

Proof. The first step is to show that t)Cj < 1/4, where u is as in f|T4l . Since 
r < 1/2, B(x°,2-f(9)r) C D' . Observe that for any x,y E D', 

\\f(x°r\f(x) - f'(y))\\ = || {fixT 1 + (fix )- 1 - fixT 1 )) 

(/'(*) - f\y)) II 

(/V)-/V))/VrW*)-/'fo))ll 

< u D ,\\x-y\\ + \\f\x*)-\f\x*)-f\x»))\\- 
\\f{x»)-\f{x)-f{y))\\ 

< lo d < \\x - y\\ + 
UJD/ \\x*-x \\-\\f(x )- 1 (f(x)-f(y))\\ 

< u) D i\\x-y\\ + 
u D/ v\\f'(x )- 1 (f'(x)-f(y))\\. (20) 

Since f lT9l implies 

l-u D ,r > 1/7(0) > 0, (21) 
the inequality (I2"0l becomes 



|/V)~ 1 (/ / (^)-/ , (2/))||< 



1 — UD>r 



\x 



I I' l l"' 

- < (22) 
1 — u D ir 

where uj is the smallest Lipschitz constant of / / ( x °) -1 /' over D 1 . 
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Recall that f(x*) = and X C D'. Observe that 

n = ||f(*T7(*°)|| 

= \\f{x»)-\f{x G )-f{x*))\\ 

< {m.a^\\f{x G )- x f\x)\\^ ■ \\x°-x*\\ 

< (max \\f'(x°)-\f'(x) - f'(x )) + f (x )' 1 f (x°)\ 

< f max ll/'Cx )-^/'^) - /'(x ))!! + lV r 
\ xex J 

< (ur + l)r. (23) 
Using pi). PH. (122]). and (1231 together give 

The last step is to to verify the other condition that B(x°, p_) C B(x°, 2j(0)r). 
Noting that y/l - 2h > 1 - 2h for < h < 1/2, it is seen that 



1- y/T^2^ 

p-ynM = 

< 2rj 

< 2(ujr+l)r 

< 2^{6)r. □ 

Next results are concerned with the size of the patch satisfying the ex- 
clusion test. 

Lemma 7.2. Let f : C n — > C n be a polynomial function with generic coef- 
ficients. Assume that all zeros of f are invertible. Let x° be a point in R n . 
U 

\\f'( x ')-if( x °)\\<J- (24) 
for all complex zeros x* of f , then there exists x* , a zero of f , such that 



■ < i^yT^jwwrmi (25) 

UJf 

= a(x*,x°). 
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Proof. By the assumption that / has generic coefficients, the polynomial / 
has a finite number of zeros. Let x*, x* 2 , . . . , x* d be all the complex zeros of /. 
Recall that a multiple zero has singular Jacobian. Hence, / has no multiple 
zeros by assumption. 

Define the polynomial f(x) = f(x) — f(x°). Note that x° is a zero of /. 
We apply the Kantorovich's theorem for complex functions (see [5J) to each x* 
with respect to /. For each x*, we use D = B(x*, p*(x*)) and u = Uf. Since 
T] = ||/'(a;*) _1 /(a;*)|| = \\f'(x*)~ 1 f(x°)\\, the assumption ([21]) guarantees 
that the condition rju; < 1/2 is satisfied. The condition B(x*, p_) C D is also 
satisfied by the definition of D. Therefore, the Kantorovich theorem states 
that there is a zero of /, call it x*, such that 

\\x*-x*\\ <<r(x*,x°). (26) 

Recall that, for any j, x* is the unique zero of / in B(x*j, p*(xj)). Therefore, 

- x* 1 1 > max{p* (x* ) , p*(x*) },i ^ j. (27) 

But (J26J) and ([27} together imply that 

x*^x*,i^j. (28) 

Hence the mapping x* x* is injective. But since / has generic coefficients 
and / and / are of the same degrees, / has at least as many zeros as / [3] . 
This implies that x° = x*, for some i. The lemma follows. □ 

Theorem 7.3. Let f(x) = f(u, v) be a polynomial system in basis <pi(u)<f)j(v) 
in two dimensions with generic coefficients. Let x° = (u°,v°) be a point in 
[0,1] 2 such thatf'( x°) is invertible and f(x°) ^ 0, x* be the closest zero in M? 
of f to x° , and 8 denote — x*\\. Let r > be such that B(x°,r) C [0, l] 2 . 
Assume 5 > Define f{u,v) such that 

,/ 2r 2/ir n 2r A 2hr n 

f{u, v) = f{j—jU ~ j—} +u° + r, —v - — + v° + r). (29) 

In other word, f is a polynomial in basis 4>i(u)(pj(v) that reparametrizes with 
[I, h] 2 the surface defined by f over the patch B(x°, r). The bounding polygon 
of {f(u,v) : I < u,v < h} satisfying itemlE of the basis properties listed in 
Section O does not contain the origin if 

r ^ 2flcond(/) 2 - (30) 
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Proof. Let X denote the patch B(x°,r) and x denote an arbitrary point in 
X. Since 5 > the contrapositive of Lemma 17.21 implies there exists a 

zero x* of / satisfying \\ f (x*)^ 1 f (x°) \\ > Therefore, the condition f[3"Ul) 
implies 

1 



r < 



< 



26 cond(/) 2 
1 



< 



More specifically, we have 



26u f \\f\x*)- l f\x) 

d\\fi{x*)-^'{x)W 



r < 



0max tf6 x||/ / (^)- 1 / , (j/)ll 



which is equivalent to 



9 ■ max /V)"7' (y) ■ r < f (*T7(^) 



y&X 



(31) 



Recall that max ye x \\f'( x *) 7'(?/)ll is the Lipschitz constant for f'(x*) 7 
on X. Hence, for any 



lArrV^-AzTV^ )!! < max||f(r)-V'(y)M|a:-a: | 



< max||f(xT7(?/)||-r. 



(32) 



Combining (13T]) and (|32|) gives 



" 7'(*T7(*) - /'(*T7(s ) < /'(sT7(* ) 



which is equivalent to 



e ■ 



fixT'm - /'(^T7>°) < /'(*T7(* U ) 



for some x G [/, /i] 2 , where x is the rescaled x and x° is the rescaled x° 
according to f l29l) . In particular, 



9 ■ max 



/V)~7(*) - /V) _ 7(*°) < fix*)- 1 f{x°) . (33) 
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Let h(x) = f'(x*)- l f(x) and g(x) = h(x) - h(x°). By §}, 

\\z\\ < 9 ■ max \\g(x)\\ , (34) 

x£[l,h] 2 

for any z in the bounding polygon P g of {g(x) : x G [Z,/i] 2 }. Since the 
bounding polygon is required to be translationally invariant (item [2] of the 
basis properties listed in Section [3]), (1341 is equivalent to 

\\y - h(x°)\\ < 9 ■ max \\h{x) - h(x°))\\ , (35) 

xG[l,h] 2 

for any y in the bounding polygon P h of {h(x) : x G [/,/i] 2 }. Substituting 
( IB!)]) into the left hand side of (13^1) yields 

\\y - h{x°)\\ < \\h{x°)\\ , (36) 

which implies that P does not contain the origin. Since fix*)" 1 is invertible 
and the bounding polygon is affmely invariant, the bounding polygon of 
{f{x) : x G [/, h} 2 } does not contain the origin, either. □ 

Theorem 7.4. Let f(x) = f(u, v) be a polynomial system in basis <pi(u)<j)j(v) 
in two dimensions with generic coefficients whose zeros are sought. Let 
X = B(x°,r) be a patch under consideration during the course of the KTS 
algorithm. The algorithm does not need to subdivide X if 

r < 1 min f 1 ~ l ±M 1 \ (3 7) 

"2 I Uf '2^cond(/) 2 J • [6{) 

Proof. If 5 > 1/wf, where S is the distance between x° and the closest zero 
x*, r < (1 — 1/7(0)) /(2uf) implies that X does not contain a zero. There- 
fore, r < 1/(4:9 cond(/) 2 ) implies that X is excluded by the exclusion test 
according to Theorem 17.31 

Observe that < uif. If 8 < — , for any x G X , 

11 *ll ^ II OllillO *|| 

\\x — x II < \\x — x I) + \\x —x II 

< r + 5 

< 1-1/7W , 1 

2 2 

< — < — = P*- 

UJf uj* 
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In other word, X is contained within B{x* l p 1f ) l a safe region and therefore 
is excluded, provided that x* is found before X is checked against all safe 
regions. By Theorem l7.ll x* is found by a region of size 2r < (1— 1/^(9)) /ud'- 
Since KTS examines larger regions before smaller ones, x* is found before X 
is checked against safe regions. □ 

It should be noted that 

1 - l/ 7 (0) > 1/(180) (38) 

hence both terms of the right hand side of fl37|) are asymptotically linear in 
1/6. The inequality fl38l) follows from the fact that 

v / TT^ < 1 + a/2 - a 2 /9 (39) 

for any a e [0, 1/4]. To prove ([39}, simplify (J39J) to a 2 - 9a + 9/4 > 0, whose 
left hand side is a convex quadratic polynomial that crosses the x-axis at 
(9 - Qy/2)/2 « .2574 and at (9 + 6v^)/2. 



8 Computational results 

The KTS algorithm is implemented in Matlab and is tested against a num- 
ber of problem instances with varying condition numbers. As Bezier sur- 
faces are widely used in geometric modeling, we choose to implement KTS 
for the Bernstein basis case. Most of the test problems are created by 
using normally distributed random numbers as the coefficients c^'s of /. 
For some of the test problems especially those with high condition number, 
some coefficients are manually entered. The degrees of the test polynomials 
are between biquadratic and biquartic. As an example, the test case with 
cond(/) = 3.5 x 10 3 is c 00 = (1.2, .5) T , c 01 = (-.6, -.6) T , c 02 = (.1,1.1) T , 
c 10 = (-1.1, -.3) T , c n = (.6,-2.3f, c 12 = (-2,-.l) T , c 20 = (.6, 1.2) T , 
c 2 i = (— 1.1, — 1.2) T , and c 22 = (— .5, A) T . This is the test problem for the 
result in the second row of Table [2] 

For the experiment, we use the algorithm by Jonsson and Vavasis [9] to 
compute the complex zeros required to estimate the condition number. Ta- 
ble [2] compares the efficiency of KTS with its condition number. The total 
number of subpatches examined by KTS during the entire computation, the 
width of the smallest patch among those examined, and the maximum num- 
ber of Newton iterations (in the cases with more than one zero) to converge 
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Num. of 


Distance 


Num. of 


Smallest 


Max. num. of 


cond(/) 


zeros 


between two 


patches 


width 


Newton 






closest zeros 


examined 




iterations 


6.0 x 10 2 


1 




21 


.0625 


3 


3.5 x 10 3 


2 


.4196 


29 


.0625 


3 


8.3 x 10 4 


2 


.6638 


33 


.0625 


3 


1.6 x 10 5 


1 




41 


.03125 


4 


2.2 x 10 7 


3 


.3624 


57 


.03125 


4 


1.3 x 10 s 


4 


.2806 


81 


.015625 


6 


1.9 x 10 9 


4 


.3069 


69 


.03125 


6 


2.0 x 10 10 


2 


.7810 


105 


.015625 


6 


2.9 x 10 11 


1 




257 


.0039 


9 



Table 2: Efficiency of KTS algorithm on problems of different condition 
numbers. 

to a zero are reported. The result shows that KTS needs to examine more 
number of patches and needs to subdivide to smaller patches as the condition 
number becomes larger. Note that the high number of Newton iterations of 
some test cases is due to roundoff error. 

9 Conclusion and future directions 

We present KTS algorithm for finding the intersections between a para- 
metric surface and a line. By using the combination of subdivision and 
Kantorovich's theorem, our algorithm can take advantage of the quadratic 
convergence of Newton's method without the problems of divergence and 
missing some intersections that commonly occur with Newton's method. We 
also show that the efficiency of KTS has an upper bound that depends solely 
on the conditioning of the problem and the representation of the surface. 
Nevertheless, there are a number of questions left unanswered by this article 
such as 

• Extensibility to piecewise polynomial surfaces and/or NURBS. 

Since KTS only requires the ability to compute the bounding polygon 
of a subpatch that satisfies the list of basis properties, it may be pos- 
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sible to extend KTS to handle these more general surfaces if bounding 
polygons having similar properties can be computed relatively quickly. 

• Tighter condition number. The condition number presented earlier 
seems overly loose. It is likely that a tighter condition number exists. 
If a tighter condition number is found, we would be able to calculate a 
tighter bound on the time complexity of KTS, too. 

• The necessity of the generic coefficients assumption. Is it pos- 
sible to analyze the efficiency of KTS without this assumption? 

• Using KTS in floating point arithmetic. In the presence of round- 
off error, we may need to make adjustments for KTS to be able to 
guarantee that all zeros are found. In addition, the accuracy of the 
computed zeros would become an important issue in this case. 

• Choice of polynomial basis. It is evident from Table CD that the 
Chebyshev basis has the best (smallest) value of 9, and therefore ought 
to require the fewest number of subdivisions. Our preliminary com- 
putational results [T7j comparing bases, however, do not indicate a 
clear-cut advantage for the Chebyshev basis. Therefore, the impact 
of the choice of basis on practical efficiency is an interesting topic for 
further research. 
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